Preparations

Load the necessary libraries

library(car)       #for regression diagnostics
library(broom)     #for tidy output
library(ggfortify) #for model diagnostics
library(sjPlot)    #for outputs
library(knitr)     #for kable
library(effects)   #for partial effects plots
library(emmeans)   #for estimating marginal means
library(MASS)      #for glm.nb
library(MuMIn)     #for AICc
library(tidyverse) #for data wrangling
library(brms)
library(tidybayes)
library(broom.mixed)
library(rstan)
library(patchwork)
library(DHARMa)

Scenario

In an honours thesis from (1992), Mullens was investigating the ways that cane toads ( Bufo marinus ) respond to conditions of hypoxia. Toads show two different kinds of breathing patterns, lung or buccal, requiring them to be treated separately in the experiment. Her aim was to expose toads to a range of O2 concentrations, and record their breathing patterns, including parameters such as the expired volume for individual breaths. It was desirable to have around 8 replicates to compare the responses of the two breathing types, and the complication is that animals are expensive, and different individuals are likely to have different O2 profiles (leading to possibly reduced power). There are two main design options for this experiment;

  • One animal per O2 treatment, 8 concentrations, 2 breathing types. With 8 replicates the experiment would require 128 animals, but that this could be analysed as a completely randomized design
  • One O2 profile per animal, so that each animal would be used 8 times and only 16 animals are required (8 lung and 8 buccal breathers)

Mullens decided to use the second option so as to reduce the number of animals required (on financial and ethical grounds). By selecting this option, she did not have a set of independent measurements for each oxygen concentration, by repeated measurements on each animal across the 8 oxygen concentrations.

Toad

Format of mullens.csv data file

BREATH TOAD O2LEVEL FREQBUC SFREQBUC
lung a 0 10.6 3.256
lung a 5 18.8 4.336
lung a 10 17.4 4.171
lung a 15 16.6 4.074
... ... ... ... ...
BREATH Categorical listing of the breathing type treatment (buccal = buccal breathing toads, lung = lung breathing toads). This is the between subjects (plots) effect and applies to the whole toads (since a single toad can only be one breathing type - either lung or buccal). Equivalent to Factor A (between plots effect) in a split-plot design
TOAD These are the subjects (equivalent to the plots in a split-plot design: Factor B). The letters in this variable represent the labels given to each individual toad.
O2LEVEL 0 through to 50 represent the the different oxygen concentrations (0% to 50%). The different oxygen concentrations are equivalent to the within plot effects in a split-plot (Factor C).
FREQBUC The frequency of buccal breathing - the response variable
SFREQBUC Square root transformed frequency of buccal breathing - the response variable

Read in the data

mullens = read_csv('../public/data/mullens.csv', trim_ws=TRUE)
## 
## ── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## cols(
##   BREATH = col_character(),
##   TOAD = col_character(),
##   O2LEVEL = col_double(),
##   FREQBUC = col_double(),
##   SFREQBUC = col_double()
## )
glimpse(mullens)
## Rows: 168
## Columns: 5
## $ BREATH   <chr> "lung", "lung", "lung", "lung", "lung", "lung", "lung", "lung…
## $ TOAD     <chr> "a", "a", "a", "a", "a", "a", "a", "a", "b", "b", "b", "b", "…
## $ O2LEVEL  <dbl> 0, 5, 10, 15, 20, 30, 40, 50, 0, 5, 10, 15, 20, 30, 40, 50, 0…
## $ FREQBUC  <dbl> 10.6, 18.8, 17.4, 16.6, 9.4, 11.4, 2.8, 4.4, 21.6, 17.4, 22.4…
## $ SFREQBUC <dbl> 3.255764, 4.335897, 4.171331, 4.074310, 3.065942, 3.376389, 1…

Exploratory data analysis

Model formula: \[ y_i \sim{} \mathcal{Pois}(\lambda_i)\\ ln(\lambda_i) =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]

where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of copper, distance and their interaction on the number of number of worms. Area of the place segment was also incorporated as an offset. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual plates.

Fit the model

Model validation

Model investigation / hypothesis testing

Further analyses

Summary figures

sqrt (gaussian) with lmer

mullens.grid = with(mullens,
   list(BREATH=levels(BREATH),
     O2LEVEL=seq(min(O2LEVEL), max(O2LEVEL), len=100)
   )
)
newdata = emmeans(mullens.lmer1a, ~O2LEVEL|BREATH,
                  at=mullens.grid, type='response') %>% as.data.frame %>%
          mutate(emmean = emmean^2,
                 lower.CL=lower.CL^2,
                 upper.CL=upper.CL^2)
head(newdata)

ggplot() +
    geom_ribbon(data=newdata,
                aes(ymin=lower.CL,ymax=upper.CL,
                    x=O2LEVEL, fill=BREATH), alpha=0.3)+    
    geom_line(data=newdata,
              aes(y=emmean, x=O2LEVEL, color=BREATH)) +
    theme_classic()

sqrt (gaussian)

mullens.grid = with(mullens,
   list(BREATH=levels(BREATH),
     O2LEVEL=seq(min(O2LEVEL), max(O2LEVEL), len=100)
   )
)
newdata = emmeans(mullens.glmmTMB1a, ~O2LEVEL|BREATH,
                  at=mullens.grid, type='response') %>% as.data.frame %>%
          mutate(emmean = emmean^2,
                 lower.CL=lower.CL^2,
                 upper.CL=upper.CL^2)
head(newdata)

ggplot() +
    geom_ribbon(data=newdata,
                aes(ymin=lower.CL,ymax=upper.CL,
                    x=O2LEVEL, fill=BREATH), alpha=0.3)+    
    geom_line(data=newdata,
              aes(y=emmean, x=O2LEVEL, color=BREATH)) +
    theme_classic()

beta

mullens.grid = with(mullens,
   list(BREATH=levels(BREATH),
     O2LEVEL=modelr::seq_range(O2LEVEL, n=1000)
   )
)
newdata = emmeans(mullens.glmmTMB2b, ~O2LEVEL|BREATH,
                  at=mullens.grid, type='response') %>% as.data.frame 
head(newdata)

ggplot() +
    geom_ribbon(data=newdata,
                aes(ymin=lower.CL,ymax=upper.CL,
                    x=O2LEVEL, fill=BREATH), alpha=0.3)+    
    geom_line(data=newdata,
              aes(y=response, x=O2LEVEL, color=BREATH)) +
    theme_classic()

mullens <- mullens %>%
    mutate(resid=resid(mullens.glmmTMB2b),
           fitted=fitted(mullens.glmmTMB2b),
           obs=fitted+resid
           )

ggplot() +
    geom_point(data=mullens, aes(y=obs, x=O2LEVEL, color=BREATH)) +
    geom_ribbon(data=newdata,
                aes(ymin=lower.CL,ymax=upper.CL,
                    x=O2LEVEL, fill=BREATH), alpha=0.3)+    
    geom_line(data=newdata,
              aes(y=response, x=O2LEVEL, color=BREATH)) +
    theme_classic()
mullens.glmmTMB = glmmTMB(pzBUC ~ BREATH*poly(O2LEVEL, 3) + (1|TOAD), data=mullens,
                          family=beta_family(link = "logit"))

#mullens.glmmTMB1 = glmmTMB(pzBUC ~ BREATH*poly(scale(O2LEVEL), 3) +
#                               (poly(scale(O2LEVEL), 3)|TOAD), data=mullens,
#                          family=beta_family(link = "logit"))

mullens.glmmTMB1 = glmmTMB(pzBUC ~ BREATH*poly(O2LEVEL, 3) + (O2LEVEL|TOAD), data=mullens,
                           family=beta_family(link = "logit"))
AIC(mullens.glmmTMB, mullens.glmmTMB1)


## library(lme4)
## mullens.glmer = glmer(pBUC ~ BREATH*poly(O2LEVEL, 3) + (1|TOAD),
##                       data=mullens, family=Beta(link='logit'))
## mullens.glmer1 = glmer(pBUC ~ BREATH*poly(O2LEVEL, 3) + (BREATH|TOAD),
##                       data=mullens, family=binomial(link='logit'))
## AIC(mullens.glmer,mullens.glmer1)
##                                         # Try this trick instead
## mullens.glmer = glmer(cbind(FREQBUC,100-FREQBUC) ~ BREATH*poly(O2LEVEL, 3) + (1|TOAD),
##                       data=mullens, family=binomial(link='logit'))

## mullens.glmer1 = glmer(cbind(FREQBUC,100-FREQBUC) ~ BREATH+poly(O2LEVEL, 3) + (1|TOAD),
##                       data=mullens, family=binomial(link='logit'))
## AICc(mullens.glmer, mullens.glmer1)

#ggplot() +
#    geom_point(data=NULL, aes(y=resid(mullens.glmmTMB), x=fitted(mullens.glmmTMB)))


#plot(mullens.glmmTMB)


plot_model(mullens.glmmTMB1, type='diag')
performance::check_model(mullens.glmmTMB1)
mullens.resid = simulateResiduals(mullens.glmmTMB1,  plot=TRUE)

plot(allEffects(mullens.glmmTMB1))
plot(allEffects(mullens.glmmTMB1),  multiline=TRUE,  ci.style='bands')
plot_model(mullens.glmmTMB1, type='eff', terms=c('O2LEVEL', 'BREATH')) 


summary(mullens.glmmTMB1)
tidy(mullens.glmmTMB1, conf.int=TRUE,  exponentiate=TRUE)

#model.matrix(~BREATH*poly(O2LEVEL,3), mullens)
emtrends(mullens.glmmTMB1,  ~BREATH,  var='O2LEVEL',  max.degree=3) 

emmeans(mullens.glmmTMB1, ~O2LEVEL|BREATH) #%>% as.data.frame

emmeans(mullens.glmmTMB1, ~O2LEVEL|BREATH,type='response')
#contrast(emmeans(mullens.glmmTMB, ~O2LEVEL|BREATH,type='response'),interaction='poly')

mullens.grid = with(mullens,
   list(BREATH=levels(BREATH),
     O2LEVEL=seq(min(O2LEVEL), max(O2LEVEL), len=100)
   )
)
newdata = emmeans(mullens.glmmTMB1, ~O2LEVEL|BREATH,
        at=mullens.grid, type='response') %>% as.data.frame
head(newdata)

ggplot() +
    geom_ribbon(data=newdata,
                aes(ymin=lower.CL,ymax=upper.CL,
                    x=O2LEVEL, fill=BREATH), alpha=0.3)+    
    geom_line(data=newdata,
              aes(y=response, x=O2LEVEL, color=BREATH)) +
    theme_classic()

r.squaredGLMM(mullens.glmmTMB1)
performance::r2_nakagawa(mullens.glmmTMB1)

References